Synchronization of coupled limit cycles 



^H Georgi S. Medvedev * 

O 
(N 

Q December 2 1 , 20 1 

Q 

Q Abstract 

< 
^ A unified approach for analyzing synchronization in coupled systems of autonomous differential 

f^ equations is presented in this work. Through a careful analysis of the variational equation of the coupled 

system we establish a sufficient condition for synchronization in terms of the geometric properties of the 

K*" local limit cycles and the coupling operator This result applies to a large class of differential equation 

C^^ models in physics and biology. The stability analysis is complemented with a discussion of numerical 

o 

■^ simulations of a compartmental model of a neuron. 

O 
O 

l; 1 Introduction 

X 

^ Synchronization of coupled oscillators has been studied extensively due to its importance in diverse prob- 

lems in science and engineering 121 [38l |42l |50l . Technological applications of synchronization include 
coordination of activity in power, sensor, and communication networks |[T4l ; and control of the groups of 
mobile agents Il39ll43i . Many experimental systems exhibit synchrony: electrical circuits [T], coupled lasers 
B4l . and Josephson junctions |[54l . to name a few. In population dynamics, the theory of synchronization 
is used to study collective dynamics Il3lll9l|50l. Synchronization plays a prominent role in physiology and 
in neurophysiology, in particular. It is important for information processing in the brain ll45l . attention. 
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arousal |[53l . and regulation of the sleep-wake cycles lITTll . In addition, synchronization underlies several 
common neurodegenerative pathologies such as epilepsy ll52]| and Parkinson's Disease ||29]| . This hst can be 
continued. 

A large body of mathematical and physical literature is devoted to different aspects of synchronization 
in differential equation models. For weakly coupled systems, very effective techniques have been developed 
Il3l]|30l[l5l[l0l[n][i6l|2l[32l|4ll (see also Chapter 10 in [27] for a survey of available methods). Synchro- 
nization in networks with symmetries has been studied in BTl (see also ||T9l and references therein). Chaotic 
synchronization has been a subject of intense research |[Tl|2l[8l|9l|2Tl|28l|40l. For strongly coupled systems, a 
number of studies explain synchronization in specific physical and biological models lini4l [T2l[T6l[37ll46ll48]| 
and for certain canonical network topologies such as nearest-neighbor coupling on a lattice ||2ll23]|. 

Networks arising in applications feature a rich variety of oscillatory mechanisms and coupling architec- 
tures. Therefore, analytical results elucidating synchronization in systems under general assumptions on the 
local dynamics and coupling operators are important. For systems with strong coupling, such results are 
rare. The exception is the work by V. Belykh, I. Belykh, and Hasler 151 |6l, where sufficient conditions for 
synchronization are given in terms of the properties of the graph of the network. The pioneering work of 
Afraimovich, Verichev, and Rabinovich 1 1 1 already identified dissipation produced by the coupling opera- 
tor as a principal ingredient in a common mechanism of synchronization. Dissipativity of the coupling was 
shown to be responsible for synchronization in many systems of coupled differential equations Iini2l l23ll46l . 
It is also relevant to stability of spatially homogeneous states in reaction-diffusion systems ||24| . However, 
the dissipation produced by the coupling alone is often not sufficient for synchronization. The analysis of 
forced Duffing oscillators coupled through position variables in f23\ shows that the intrinsic dissipativity of 
the individual (local) subsystems is equally important. In fact, synchronization is achieved by the interplay 
of the dissipativity of the coupling and the intrinsic properties of the local systems. The analysis in i231 
describes an important mechanism of synchronization using the Lyapunov functions that are specific to the 
Duffing systems coupled via a discrete Laplacian. The goal of this paper is to study this mechanism under 
general assumptions on the local oscillatory dynamics and for a broad class of coupling schemes. Through 
a careful analysis of the variational equation of the coupled system, we derive a sufficient condition for 



synchronization in terms of the geometric properties of the local limit cycles and the coupling operator. To 
achieve this, in the vicinity of the periodic solution of the coupled system, we construct a moving frame 
of reference. After a series of coordinate transformations, we arrive at a system of equations that reveals 
the interplay of coupling and local dynamics, and shows their combined contribution to the stability of the 
synchronous solution of the coupled system. 

The key step in this analysis is finding a suitable transformation for the coupling operator in the moving 
coordinates. As a by-product, we develop a rigorous reduction of the coupled system to the system of 
equations for the phase variables. In approximate form, this system of phase equations was obtained in [34]. 
Similar to Kuramoto's phase reduction for weakly coupled systems, we expect that these phase equations 
will be useful in studies of physical and biological coupled oscillator models in the strong coupling regime. 
For analytical convenience and because the coupled limit cycle oscillators are common in applications, in 
this paper we consider local systems whose dynamics are generated by limit cycles. Synchronization of 
chaotic systems as considered in ITJ |2l |23l is admittedly more appealing and physically less intuitive than 
synchronization of coupled limit cycles. However, from an analytical point of view, the latter problem 
contains many of the ingredients responsible for synchronization of systems with more complex dynamics. 
For a discussion of how the results of the present study can be extended to chaotic synchronization and for 
related extensions for systems with time-dependent and nonlinear coupling schemes and randomly perturbed 
local systems, we refer the interested reader to ll34ll . 

The paper proceeds as follows. In Section[2]we list our assumptions and state the main result. Section[3] 
explains the assumptions made in the previous section in the context of three examples. This is followed by 
the proof of the main theorem in Section |4] 



2 Assumptions and the main result 

2.1 The local dynamics 

We start by discussing the assumptions on a local system: 

x = f(x), x:M^M" (2.1) 

where f : R" — )• M" is continuous together with partial derivatives up to second order. We assume that 



X = u(t) is a periodic solution of (2.1 1 of period 1 with a nonvanishing time derivative u(t) 7^ 0, t € S-*^. 
Denote the corresponding periodic orbit O = {x = u(t), t G S^ := M^/Z}. Near O, one can introduce an 
orthonormal moving coordinate frame (cf. Theorem VI. 1.1, ||22]| ): 

{v(e),zi(0),Z2(0),...,Zn-i(0)}, v(^) = ^, eGS\ (2.2) 

The change of variables 

x = u{e) + z{e)p, z{e) = coi(zi(^), . . . ,zn_i(e)) e m"x("-i) (2.3) 

in a sufficiently small neighborhood of O, defines a smooth transformation x i— ;• {9, p) £ S^ x M"~^ ||22l . 



Lemma 2.1. In new coordinates (2.3 ), near O local system (2.1 ) has the following form 



e = i + z'{e)p + o{\p\^), (2.4) 

p = k{e)p + 0{\p\^), (2.5) 



where 



A(^) = z'{e)D^{u{e))z{9)-z'{e)z'{e), {i.i) 



where M^ stands for symmetric part of matrix M, M^ := 2 "^(M + M ). 

Notational convention. To simplify notation, in the calculations through out this paper, we will often 
suppress 9 and u{9) as arguments o/f, Z,v, etc, when the expression of the argument is clear from the 
context. We continue to denote the differentiation with respect to t and 9 by dot and prime respectively. 



Proof. The proof follows the lines of the proof of Theorem VI.1.1 [225. By plugging (2.3 1 into (2.1 1, we 
have 

[u'{9) + Z'{9)p] 9 + Z{9)p = f (u(0)) + D^{u{9))Z{9)p + Qi(0, p), (2.8) 



where 



Qi(0, p) = f(u(0) + Z{9)p) - f(u(0)) - D^{u{9))Z{9)p = 0{\p\ 



By multiplying both sides of (2.8 1 by v |f | , we obtain 



1,,T/ 



^(1 + |f|-^' ZV) = 1 + |f|-\' (Df)Zp + 0(|p| 



(2.9) 



For small |/9|, (2.9) can be rewritten as 



^ = l + |frV((Df)Z-Z> + 0(|p|2). 



(2.10) 



By differentiating both sides of \i^ {9)Z{9) = 0, we have 



T/nJ:^T- 



v'Z' = -(v')'Z = -v'(Df)'Z. 



(2.11) 



By plugging in (2. 11 1 into (2. 10 1, we arrive at (2.6 1. Next, we multiply (2.8 1 by ZJ , and using ^ = 1 + 0(|p|) 



we derive \1.5) 

D 



Next we formulate our assumption on the stability of the local limit cycle. 



Assumption 2.2. Let fii{9) denote the largest eigenvalue of /K^{9) (cf. (2.7)). We assume 



fii{u)du =: — /i < 0. 



(2.12) 



Assumption |2.2| implies exponential stability of the trivial solution of the linear periodic system 



P = A(t)/9. 



(2.13) 



For convenience of the future reference, we formulate this statement as a lemma. 



Lemma 2.3. Suppose A(t) is a continuous periodic matrix of period 1. Then (2.12) implies that for any 
<e < n, 

\R{t)R-^{s)\ < Ci exp{-{fi -e){t- s)}, t > s, (2.14) 



where Ci > and R{t) stands for the principal matrix solution of (2.13 ) 



Proof. Let 



By periodicity of ^1(6') and (2.121, 



1 r 

u{t) ■■= - I ni{s)ds. 
t Jo 



lim u{t) = -// < 0. 

t— >CXD 



(2.15) 



(2.16) 



Further, 



d^ 

dt 



\p{t)\' = 2p'A%t)p<2f,,{tMt)\^ 



and, by Gronwall's inequality. 



\p{t)\<\pmexpMt)t}. 



(2.17) 



This in turn implies that for large t » 1 and arbitrary < e < —p 



mt)pm 
\pm 



<ex.p{-{p-e)t}. 



(2.18) 



On the other hand, by the Floquet theorem, 

R{t) = P{t)ex.-p{tA}, (2.19) 



where P{t) and A are periodic and constant matrices respectively. By comparing (2.181 and (2.191, we 



conclude that all eigenvalues of A must have negative real parts. This yields ( 2. 14 1 

D 

2.2 The coupling operator 



By the coupled system, we call a collection of N local dynamical systems (2.1 1 interacting via a linear 
coupling operator D : M^" -^ M^" 



x = f{x) + gDx, x= [x'^') ,x('> , . . . ,x('^> j G M^'", (2.20) 

where f{x) = (f(x(-^)), f(x(2)), . . . , f (x^'^^)) G M^" and <? > is a parameter controlhng the strength of 
interactions. In this paper, we consider separable schemes ||34l: 

L» = D (» L, (2.21) 

where L G M"^", D G M^^^, and ig) stands for the Kronecker product f25\i . From the modeling point of 
view, separable coupling is natural as it reflects two levels of the network organization. The global network 
architecture of interconnections between local systems is captured by D. Matrix L reflects the organization 
of the coupling at the level of a local system: it shows what combination of local variables participates 
in the coupling. The separable structure of the coupling translates naturally into the stability analysis of 
the synchronous solution, revealing what features of the global and local organization of the coupling are 
important for stability. 



Definition 2.4. By synchronous periodic solution of {2.20) (when it exists) we call 



2; = n(t) := In ® u(t), In := (1, 1, • • • , 1)' G M'\ (2.22) 



where u(t) is a periodic solution of the local systems (2.1). 



Next we specify the structure of the separable couphng operator D (cf. (2.21 1). There are two (sets of) 
assumptions. The first simpler assumption ensures that the coupled system admits a synchronous solution. 
The second set of assumptions guarantees the stability. As far as existence is concerned, we need to postulate 
that In G ker(D). We further assume that ker(D) is one-dimensional to limit our study to connected 
networks. Thus, we assume 



D G /C = {M G M^^^ : ker M = Span {In}} 



(2.23) 



Next, we turn to assumptions that ensure stability of the synchronous solution. To this end, we define an 
{N -1) X N matrix 



/ -1 1 
-1 1 








\ 



(2.24) 



y ... -1 1 y 

In the stability analysis of the synchronous solution, we will use matrix D G m(^^i)^(^^i) defined by the 
following relation 

SD = DS (2.25) 



For any D G /C, D is well-defined (cf. ||33]| . see also Appendix in II35II ). The spectrum of D is important for 
synchronization. This motivates the following definition. 



Definition 2.5. Matrices from 



2? = |m G /C : x'^Mx < Vx G M^"V{0}} 



(2.26) 



are called dissipative. 



Finally, we state our assumptions on matrix L describing the local organization of the coupling . 



Assumption 2.6. L is (symmetric) positive semidefinite, i.e., L = L and x Lx > Vx G R". 



The following definition distinguishes two important cases that come up in the stability analysis of the 



synchronous solutions of (2.20 1 and (2.21 1. 



Definition 2.7. IfL is positive definite, we say that the coupling is full (rank), otherwise we call it partial 
(rank). 



Remark 2.8. In a slightly different form, the full and the partial coupling schemes were introduced in II23II . 



Dissipative matrices yield synchronization when the coupling is full and sufficiently strong. Synchro- 
nization in the partially coupled systems requires an additional hypothesis. 

Assumption 2.9. T/'dimker L = I > 0, let {pi, p2, . . . , Pi} be an orthonormal basis o/ker L Denote 
orthogonal matrix = col(v, zi, Z2, . . . , Zn-i) and matrix 



L!(t) 



\c{t) vT(t)(Df(u(t))rZ(t) 

ZT(t)(Df(u(t)))=v(t) A^(t) - lc(t)l._i 



where c = v"'"(Df)v. Define matrix 



(2.27) 



G(t) = (gij(t)) G R'><^ gij(t) = (Ll(t)0"r(t)pi,0"r(t)pj). 



(2.28) 



Let Ai(t) denote the largest eigenvalue of symmetric matrix G(t). We assume that 



I Ai(t)dt<0. 
Jo 



(2.29) 



Assumption 2.9 captures the geometric properties of the local limit cycle through matrix L^ (cf. (2.27 1), 



which depends on the basis functions {v, zi, Z2, . . . , Zn_i} defined in the vicinity of the limit cycle and the 
Jacobian Df (u(t)) evaluated along the periodic orbit and the properties of the coupling operator through the 



basis for the kernel of L, {pi, p2, . . . , pi}, entering the definition of G (cf. (2.28 1). Condition (2.29 1, there- 
fore, reflects the interplay of the intrinsic properties of the limit cycle and the coupling operator and their 



contribution to the stability of the periodic solution of the coupled system. Having reviewed the assumptions, 
we state the main result of this paper. 



Theorem 2.10. Suppose a periodic solution of local system {2.1 1 is stable in the sense of Assumption 2.2 (cf 



(2.12)) and the coupling operator in (2.20) is separable (cf. (2.21)) and such that D G D and L is positive 
semidefinite. If the coupling is partial, in addition, we assume that ^2.29) holds. Then for sufficiently large 



g > 0, the coupled system (2.20) has an exponentially stable synchronous limit cycle. 



3 Discussion and examples 



We precede the proof of Theorem 2.10 with a discussion of the assumptions and the implications of the 
theorem. We illustrate our techniques with three examples. In the first example, we use a planar vector 



for which Assumption 2.9 can be checked by a simple explicit calculation. The second example is used to 
explain the assumption that the coupling operator D is dissipative and to elucidate the range of networks 
covered by this assumption. The third example is meant to illuminate the distinction between the full and 
partial coupling schemes in the context of a biophysical model of a neuron, and to show how one verifies 



the assumptions of the Theorem 2.10 in practice 



3.1 Coupled radially symmetric oscillators 



Consider radially symmetric local vector field: 



x = f(x), f(x) 



Xi -X2 -Xi(xf +xl) 
Xl +X2 -Xl(xf +X2) 



X= (X1,X2)"^ G 



(3.1) 



For coupled local systems (3.1 1 Theorem 2.10 yields the following sufficient condition for synchronization. 



Corollary 3.1. Coupled system ( 2.20\ with separable coupling (2.21 ) and local vector field (3.1 ) has an 
exponentially stable synchronous limit cycle if D & T> and nonzero matrix L is positive semidefinite. 
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Proof. If L is positive definite then Corrolary 3.1 follows from Theorem 2.10 Suppose that L has a ID 
kernel spanned by 

P= (Pl,P2)"^- 



Let 



T 



u(t) = (cost, sin t) , v(t) = (— sin t, cost) , z(t) = (—cost, — sin t 



Then 







cos t — sin t 



sin t cos t 



and L 




-2 



Further, q = O^p and 



Ai(t) = q(t)'L'iq(t) = -2(pi cost + p2 sin t)^ 



and 



Ai(t)dt = -(p^ + p^)<0. 



D 



Example 3.2. Let f = (fi, f2)^^ as in {3.1) and consider 



Ai) 



^i) 



N 



h{xf,xf), i = l,2,...,N. 



In this example, 



1 







is positive semidefinite. Suppose hij = hji > and denote 



D = (dij), dij=<^ 






Tjfdim ker D = 1 then Y) £ D, by Gershgorin 's Theorem. In general, D G P does not have to be symmetric. 
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For a more complete description of dissipative matrix we refer the reader to / l?3l/ (see also Theorem 3.4 
below). 

3.2 The consensus protocol 

To elucidate what features of the coupUng are important for synchronization we choose consensus protocols, 
a framework used for modeling coordination in the groups of dynamic agents [43 1. The continuous time 
variant for this problem deals with N agents whose states are given by x^'^\t), i = 1,2, ... ,N and are 
governed by the following system of differential equations: 

N 

x» =^dij(x(-'')-x«), i = l,2,...,N. (3.2) 

Here, weights dij ,i i^ j (which for simplicity we take constant) describe the interactions between two 



distinct agents x^*) and x'^^\ After setting da = — J2jj=i ^ij' ^^ rewrite (3.2| in a vector form 



± = Dx, D = {dij) elC, x = (x(^\2;(2),...,x(^)). (3.3) 

The problem data can be conveniently represented by a weighted graph Q = (V, £,{dij}), where vertex set 
V = {1,2,..., N} lists all agents, the pairs of interacting agents are recorded in the edge set £, and the 
weights {dij} quantify the intensity of interactions. If Q is connected ker D = SpanjlNJ. Both positive 
and negative weights {dij} (corresponding to synergistic and antagonistic interactions respectively) are 
admissible. Likewise, the interactions do not have to be symmetric, i.e., dij may differ from dji. In designing 
consensus protocols, one is interested in weighted graphs yielding convergence to spatially homogeneous 
state 

|xW(t)-x(^')(t)| ^Oast^oo. (3.4) 



The following questions related to ( |3.3| ) are important in applications: determining the rate of convergence 
in (|3.4[) and relating it to the network topology and weight distribution Il55ll56l . finding an optimal weight 



distribution yielding (3.3 1 a desired property such as the fastest convergence (under certain constraints) 
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Figure 1: Schematic representation of spatial structure of a compartmental model: (a) linear cable, (b) 
branched cable. Dynamic variables v^*) and u^'^\ i = 1,2, . . . ,N, approximate voltage and calcium concen- 
tration in each compartment. 



55]|5Tj or minimal effective resistance 1201 . or robustness to noise 11561 . to name a few. For studying these 



questions it is important to describe all weighted graphs that endow (3.4 1 with synchronous dynamics. For 



the discrete time counterpart of (3.2 1, this question was answered in ||55]| . The following theorem describes 



all such graphs for the continuous time problem ( 3.2 1. 



Theorem 3.3. Let x{t) be a solution of the initial value problem (3.3) with initial condition xq G M . Then 
dl^l) holds for any xq^R^ iffD e V. 



Proof. By multiplying both sides of (|3.3|) by S (cf. (2.24 1), we have 



y = Dy, y := Sx. 



(3.5) 



The equilibrium of (3.5 1 is asymptotically stable iff the symmetric part of D is negative definite, i.e., when 
D eV. 

D 



Therefore, dissipative coupling matrices are precisely those that enforce synchrony in ( |3.3[ ). Remarkably, 
dissipative matrices admit an explicit characterization. 



Theorem 3.4. [34^] T> eViff 



D = QAo 



(3.6) 
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for some Q G 



dTVxTV 



with negative definite symmetric part and 



Ao := S^S 



/ 1 -1 ... 
-1 2 -1 ... 







-1 1 



(3.7) 



Theorem 3.4 gives a convenient computational formula for D: 



D = SQS^. 



(3.8) 



This formula can be used for studying the rate of convergence of solutions of ( 3.2 ) to the homogeneous state 



for different network topologies. For some canonical network architectures, including nearest neighbor and 
all-to-all coupling schemes, the rate of convergence (i.e., the largest eigenvalue of D) can be found analyt- 
ically for other, such as networks with random connection weights, the rates can be computed numerically 
using p.8[). We refer an interested reader to |[34ll33l . where these examples are discussed in detail. 



Theorem 2.10 extends the argument used in Theorem 3.3 to networks whose local systems have multidi- 



mensional phase space and nontrivial dynamics. In the analysis of the general case, in addition to the global 
network architecture (i.e. D G V), the coupling organization on the level of the local systems becomes im- 
portant. For instance, it matters what local variables and in what manner are engaged in the coupling. This 
information is contained in matrix L (cf. ( 2.21| )). The analysis of the separable coupling schemes highlights 
the importance of the full versus partial coupling distinction. In the latter case, synchronization depends on 



the combination of the properties of the local dynamics and the coupling operator (cf. Assumption 2.9 1. The 



example in the following subsection is chosen to provide a reader with an intuition for the mechanisms of 
synchronization in fully and partially coupled systems. 
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3.3 The compartmental model 



The following systems of differential equations is a nondimensional model of the dopamine neuron (cf. 

my. 



ev 



(i) 



<7i(7;«) (i?i -t;«) +g2{n^''>) [e2-v'^'^) +53 [e, - v^''') + 1^\ 



;,i^ 



u 



uigi{v'^'^)(Ei-v 



,W 



,(0 



+ 



/», i = l,2,...,Ar. 



(3.9) 



(3.10) 



Here, v^'^' and u^''> approximate membrane potential and calcium concentration in Compartment i of an 
axon or a dendrite of a neuron. Equations ( |3.9| ) and ( 3.10| ) describe the dynamics in each compartment 
using Hodgkin-Huxley formalism (see ||T3]| for more background on compartmental models). The nonlinear 



functions gi{v) and 52 (^) and the values of the parameters appearing on the right hand sides of (3.9 1 and 



( 3.10| ) are given in the appendix to this paper. Equations (3.9l and (3.10l reflect the contribution of the 
principal ionic currents (calcium and calcium dependent potassium currents) to the dynamics of v^''' and 
^t(*^ In addition, the coupling terms 



I« = flVd,.(t>(^^)-?;»^ 



N 

j=0 

N 



I« = 55^d,,(u(^- 



ij) _Jih 



(3.11) 
(3.12) 



j=0 



model the electrical current and calcium diffusion between the adjacent compartments. The off-diagonal 
entry gdi^ of matrix gD corresponds to the conductance between Compartments i and j. The structure of 
the coupling matrix D, i.e., the pattern in which nonzero entries appear in D, reflects the geometry of the 
neuron. In the simplest case of a uniform linear cable with no-flux boundary conditions (see Fig. [T^), the 
coupling matrix D = — Aq (cf. ( |3.7| )). D may have a more interesting structure, e.g., in the models of 
dendrites with more complex spatial geometry (see Fig.[Tj5). As follows from (3.9 1 and (3. 10 1 the coupling 

is separable with 

/ 
1 

61 



(3.13) 



V 
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Figure 2: Numerical simulations of compartmental model (3.9l-( 3?T0| ). (a) Phase plane plot of the limit 
cycle of the local system, trajectories of five uncoupled oscillators are plotted for (b) Timeseries of five 
coupled oscillators, (c) The largest eigenvalue /x(t) (dashed line) and integral Jq \i{ff)dQ are plotted over one 
period of oscillations. The eigenvalue takes both positive and negative values; but since the integral over one 
period is negative, the local limit cycle is exponentially stable. (d)The largest eigenvalue Ai(t) (dashed line) 



of matrix G(f ) determining stability of the synchronous regime when the coupling is partial (cf. (2.28 1). In 
solid line we plot Jg Ai (t)dt for one period of oscillations. The integral of Ai (t) over one period is negative. 
Therefore, the synchronous solution is stable. 



where b\= g (5. If 5 > the coupling is full. If calcium diffusion is ignored (5 = 0) the coupling becomes 



partial. Below, we discuss the assumptions of Theorem 2.10 in relation to the model at hand. 



We start with the conditions on the coupling matrix D. To be specific, we assume the nearest-neighbor 
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coupling (see Fig.[T^), i.e., D = — Aq. Then 



/ 



D 



-2 10 
1 -2 1 








V 







,{N^l)x(N-l) 



J 



Clearly, D G D because D is (symmetric) negative definite. The conditions on L for both full and partial 

coupling cases obviously hold. As is typical for conductance-based models of neurons, away from the 

singular limit e — )• 0, the analytical estimates on the eigenvalues of the variational equations such as in 

Assumptions 2.2 and |2.9[ may be difficult to derive. We verify ( |2.12 ) and ( 2.29| ) numerically, after we 

briefly review basic numerics for (3.9) and (3. 10 1. Fig. [2^ shows the limit cycle of a single uncoupled 

oscillator in (3.9 1 and ( 3.10| ). The corresponding time series are shown in Fig.[2V). In Fig.[2b, we plot /Ui(0), 

the largest eigenvalue of A^(0). In this example A{6) is scalar. Note that over one cycle of oscillations 

fj-i{9) takes both positive and negative values. However, since the integral over one cycle of oscillations 
fperiod 



/o fii{6)d6 is negative (Fig. 2:), the limit cycle is exponentially stable (cf. Lemma 



2.3 1. Therefore, 



for the fully coupled variant of (3.9 1 and (3.101 all conditions of Theorem 2.10 hold. Note that in the full 
coupling case, synchronization is determined solely by the properties of the coupling operator. The only 
information about the local system which we use is the exponential stability of the limit cycle. Conditions 
for synchronization in partially coupled systems use the information about the local dynamics (through 



matrix Li (cf. ( 2.27 1)) and about the coupling operator (through ker L). The numerical verification of ( 2.29 1 
for partially coupled variant of ( 3.9 1 and ( 3. 10 1 is given in Fig. [2|i. 



4 The proof of Theorem 2.10 



The proof proceeds as follows. In ^4.1[ we construct a suitable system of coordinates near the periodic orbit 



of the coupled system. In ^ 4.2 we analyze the linear part of the variational equation. In ^ 4.3 we extend the 
stability analysis to the full system. 
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4.1 The local coordinates for the coupled system near the limit cycle 



The moving coordinates for the coupled system are obtained by combining the coordinates (2.3 1 for local 
systems 



iNn 



( v(l) \ 



.(2) 



9 X 






i^ x(N) y \ u(0(N)) + z(eW)pW y 



(0,p)g(S^)^xM(^-i)". ^4_j) 



where 6* = (6'(^) , 6'(2) , . . . 6i(^') ) , p = (/^(i) , p(2) , . . . , p(N) ) . By following the steps of the proof of Lemma|2j] 
for each local system we have 



^« = 1 + 3^(0(0)^(0 + ;^^CT + 0(|p«|2) 



|f(e(0)| 



p(0 = A(0(O)p(O+zT(0(O)CT + O(|/ 



p(0|2), i = l,2,...,N, 



(4.2) 
(4.3) 



where CT stands for the coupling terms 



N N 

CT := y^dijL (u(0(O) - u(0(J))) +5Edij'- (^(^^'V^'^ " Z(0«)p« 



(4.4) 



In the moving coordinate frame, the coupling operator becomes nonlinear. We linearize it using the Taylor's 
formula. It will be convenient to have Taylor's coefficients appearing in the expansions for local systems 
evaluated at a certain common value. For this purpose, we use the average phase defined by: 



N 



^ ■= ^Tq ^here i G ker D"^, ^ ^i = 1. 



(4.5) 



1=1 



The existence of ^ G R^ with the properties specified in (4.5 1 follows from the following considerations. 
First, because rank D^ = A^ — 1, there is nonzero ^ G kerD^. ^ can be chosen to satisfy the second 
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condition in (4.5 1 provided ^^1n / 0. Supose the contrary. Then, 



^ G (1n)^ = (kerD)^ = i?(D^) ^ 3 ?? / : C = D^ry. 



Note that ^ and r] are hnearly independent. Since ^ G ker D^, 



(D^ yr] = =^ rank (D^) <N -2. 



Therefore, either the geometric multiplicity of the zero eigenvalue of D is greater or equal to 2, or the size 
of the block corresponding zero eigenvalue in the Jordan normal form of D is greater or equal to 2. Either 
statement contradicts the assumption that D eV (see Lemma 2.5 Il33l ). 



Similarly, we define 



Q := {^ ^ \„-lf p. 



(4.6) 



From the definitions of -d and g, we have 



N N N 

0(0-^ = ^e/')-^^/J) = ^ej(0«-0«) = 0(101), 
j=i j=i j=i 

N 

j=i 



(4.7) 
(4.8) 



where cj) = SO, r = (S (8) ln)p, and S is defined in (2.24i 



Using (4.7 1, we represent 



u(0(j)) = u(T?) + u'(^)(^(J)-^) + O(|0|2), 
Z(0(j)) = Z(i9) + Z'(t?)(#-^) + 0(|(/.|2). 



(4.9) 
(4.10) 



By plugging ( 4.9 1 and ( 4. 10 1 into ( 4.4 ) and using ( 4.7 1, we have 



N N 

CT = 55;dijLu'(^)(0«-0«)+5EdijLZ(^)(p«-p(J))+o(5|0|2^5|^||^|). (4.11) 

j=i i=i 
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Using (|4.11|), we rewrite dOl and dO) 



N 



N 



ei') 



^(0 



1^1 u 

N N 

i=i j=i 



1 + aT(^)^W + 5 ^ dijm(^)0« + fj^Yl dijUiJ Wp^^ + 0(|p|', 9\<^\^9\pM), (4-12) 

(4.13) 



where 



m(??) = v^(t9)Lv(??), U(t?) = v^(i?)LZ(t?), and M(??) = Z^(t?)LZ(i9). 



(4.14) 



The averaged system is obtained by multiplying equations for O^'"' and p^^' by ^i and adding them up: 



^ = l + aT(^)^ + 0(|p|2,5|c/>|2,5| 

g = A{^)g + 0{\p\\g\(^\^g\pM,g\Q\\r\). 



(4.15) 
(4.16) 



Using tD as an independent variable, we rewrite (4.121 and (4.13 1 as follows 



N 



N 



i=i '•' ' i=i 

+ 0(|pp,5l'/'P,5|p||</'|,ffl^l|r|), 



dp« 



AT 



TV 



AWp«+5l/lEdijUj(^)^^^^+5EdijMij(^)/'' 



(J) 



+ 0{\p\',g\(t>\',g\pM,9\r\\Q\), i = l,2,...,N, 



(4.17) 



(4.18) 



or in matrix form for x := {e^^\ p^^\ 0^'^\ p'^^\ . . . , e^^\ p^^'^) 



dx 



In «> I I + 

On-1 

|2 „|i|2 



In ' 



a' 



O A 



+ gD. 



m |f|-iU 
IflU"^ M 



+ O{\p\',g\^\\g\p\\4>\,g\r\\0\), 



(4.19) 



where On-i := (0, 0, ... , 0)^ G M^"^. By multiplying both sides (4.19l by S (g) !„ and recalling (p = S9 
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andr = (S O In-i)/?, we obtain the system for y := {(l)^^\ r'^^\ (l)^^\ r^^\ . . . , (t)'^^''^\ r^^'^'^] 



dy_ 



In 



a^ \ . / m |f|-iU 

\+gT>®\ 

O A J \ |f|UT M 

+ Oi\p\^g\cP\^g\p\\cPlg\r\\g\). 



(4.20) 



The final coordinate transformation in this series is used to make the matrix multiplied by D in (4.20 1 
symmetric: 

^ = V\T\<P and z = ^=. (4.21) 

Note, 



M 
dz 






(4.22) 
(4.23) 



Using (4.21 1, (4.221, and (4.23 1 from (4.20 1 we obtain 



dy 
dd 



B{^)y + Qiiy, g), B{^) -.= gBo{^) + Bi{§), 



(4.24) 



where y:=(^W,z(i),V('),z(2),...,^(N-i),z(N-i)),Qi(y,^) = 0(|e|2,<7|y|',ff|y||e|) and 



Sn = D 



m U 
U^ M 



Bi = In-1 



I 



\c 2vT(Df)Z 

A-idn-l 



(4.25) 



and c(i?) = v"''('!9)Df(u('!9))v(7?). We complement (4.21 1 by the equation for q 



dg 



A{'&)g + Q2{y,Q), Q2{y,p) = 0{\g\\g\y\^,g\y\\g\). 



(4.26) 
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4.2 The linearized system 



In this subsection, we study the linearization of (4.24 1 



y = B{t)y. 



{A.ll) 



In the following lemma, we prove that y = is exponentially stable solution of (4.27 1. 



Lemma 4.1. Let <I>(t) denote a fundamental matrix solution of (4.27) and $(t, s) := <I>(t)<I> (s). Then 



|$(t, s)\<C2 exp{-A(t - s)}, t > s 



(4.28) 



for certain positive constants C2 and A. 



The proof of Lemma 4. 1 follows from Lemmas 4.2 and 4.3 which we prove first 



Lemma 4.2. The spectrum of symmetric matrix 



U{0) 



m{9) [J{9) 
UT(^) M{9) 



vT(0)Lv(0) y'^{6)lZ{9) 
Z"r(6')Lv(0) ZT(^)LZ(6') 



(4.29) 



coincides with that of L. In particular, for every 9 £ S^, M (9) is a positive semidefinite matrix whose rank 
is equal to rank(L). 



Proof. This follows from 



M{9) = 0'^{9)lO{9), 



where = col (v, zi, Z2, . . . , Zn-i) is an orthogonal matrix. 

D 

Lemma 4.3. Let Ai(t) denote the largest eigenvalue of 



B%t)=gB^,it) + Bm. 



(4.30) 
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Then there exists go > and A > such that 



I Xi{t)dt 
Jo 



-A, g > go- 



(4.31) 



Proof. Matrix Bq (t) = D® (g) M (t) is negative semidefinite for all t > 0, because tiie eigenvalues of D^ 
are negative and tliose of M are nonnegative. Denote the (constant) eigenvalues of Bq by A^ : 



o>a(°)>a(°)>..->a 



(0) 
{Af-l)n- 



For small 6 > 0, the eigenvalues of Bq + 6B( perturb smoothly US 



XkAt) = 4°^ + ^^k\t) + 0(5^), k = l,2,...{N- l)n. 



(4.32) 



We consider the full rank coupling case first. If L is full rank then so is M. Denote 



A := -O.bXf^ < 0. 



Choose (^n > such that 



max max A^ sU) < —A for < (5 < Sq. 

k te[o,i] 



Then ior g > gQ := 6q , the eigenvalues of B"^ = gB'Q + Bf are negative and are bounded from zero by 



gX. This shows (4.31 1 for the full coupling case. 



It remains to analyze the case ker L = / > 0, i.e., 



Al' 



(0) 



4Vi) = 0, and A^^ViH^ =: 2A < 0. 



By (4.32), 



XkAt) = ^^k\t) + 0(.S^), k = l,2,...{N- 1)1. 



(4.33) 
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For small 5 > 0, we have 



1 r-l 

max Afc5(t)(it<(5 / max Xksit)dt + 0(5 

ke{l,2,...,(N-l)n} ' Jq ke{l,2,...,{N-l)l} ' 



Thus, to show (4.31 1 we need to verify 



/ max Afc si't)dt < 0. 

Jo ke{l,2,...,{N-l)l} ' 



(4.34) 



For this, we review the construction of the correction terms A^ (t) (cf. Appendix lITSl ). Choose an or- 



(1)/ 



thonormal basis for keri?Q(i) {r/i(t),r/2(t), . • ■'ni{N-i)ii)}- Then A^ (t) are the eigenvalues of G = {g. 



»j/' 



gij = {Bl{t)rii{t),rij{t)). Below we show that Assumption 2.9 guarantees (4.34i. To this end, we con- 
struct a basis for ker i?f (t). Recall that {pi, p2, . . . , pi} stands for the orthonormal basis of ker L and 
= col(v, zi, . . . , Zn-i). We choose 



?7i(t) = ei, ®ei2(t), i = ii2-l)l + k, h G {1,2, . . . ,iV - 1}, 12 £ {1,2, . . . ,1}, 



(4.35) 



where ^i(t) = 0"'"(t)pi, ej = (5j, . . . 7<5(7v-i))^ ^ ^^ ^' ^"^^ ^i denotes the Kronecker delta. Vectors in 



(4.35 1 form an orthonormal basis of ker -Bf (t). Further, 



Qijit) = {Blit)r]iit),r]j{t)) = ((In-1 ^)(ei, ® 4(t)),ej, ® ej2(t)) = (ei, <E) 11^^,(1) , e^, (g) ^,,(1)) 
(LiCi,(t),ej2(t)) = (LiOT(t)pi„OT(t)pj,), ii =ji, 
0, otherwise. 



Thus, G = In-1 ® G where G is defined in ( |2.28 ). The eigenvalues of G are those of G taken with 



multiplicity (A^ — 1). Clearly, Assumption 2.9 is equivalent to (4.34 1. This concludes the proof of the 
lemma. 

D 



Proof. (Lemma 4.1 1 The statement of Lemma 4.1 follows from Lemma 4.3 and Lemma 2.3 



D 
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4.3 The endgame 



To complete the proof of stability of the synchronous solution, we study the initial value problem for (4.24 1 
and ( |4.26| ) 



X = d\ag{B{t),A{t))x + Q{x), 
x{0) = xo, 



(4.36) 
(4.37) 



where by abusing notation we denote x := {y, g)'^ and the independent variable by t. Matrices A(i) and B{t) 
are defined in {2.1) and (4.24i. The nonlinear terms are collected in Q{x) := {Qi{x), Q2{x)) = 0{g\x\'^). 



We suppress the dependence of Q on g, because once it is chosen sufficiently large g will be considered 
fixed. Let X{t) denote a principal matrix solution of the homogeneous system 

X = d\ag{B{t),A{t))x. 



By Lemmas 2.3 and 4.1 for X{t,s) = X{t)X-^{s) wehave 



\X{t, s)| < Cs exp{-K{t - s)}, t- s>0, 



(4.38) 



for some C3 > and < k < min{/i, A} (cf. (2.14i and (4.281). 



Lemma 4.4. Let g > go be fixed. Then for any sufficiently small e > (possibly depending on g) and any 
initial data 

\xo\<e (4.39) 



solution of the initial value problem (4.36) and (4.37) satisfies 



sup \x{t)\ < eexp < 1 > . 

t>o L 2 J 



(4.40) 



As in the famous theorem of Lyapunov on stability by the linear approximation, the statement of 



Lemma 4.4 follows from the stability of the linear part of (4.36 1 captured by (4.38 1. The proof of the 
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lemma relies on a weaker statement, which we prove first. 



Lemma 4.5. Under the assumptions of Lemma 4.4 we have 



sup \x{t)\ < 2e. 



(4.41) 



Proof. Fix 5 such that 



< 5 < min 



1 K 
2^' 2 



Because 



for sufficiently small e > we have 



Q(0) = and ^ = 0, 
ox 



(4.42) 



(4.43) 



\Q{x)\<5\x\, |j;|<2e, 
\Q{x2) - Q{xi)\ < 5\x2 - xi\, |xi_2| < 2e. 



(4.44) 
(4.45) 



Consider a functional sequence defined by 



Xl{t) = Xo, 

Xn+i{t) = X{t)xQ+ j X{t,s)Q{xn{s))ds, n = 1,2, 

Jo 



(4.46) 

(4.47) 



We use induction to show that 



sup |x„(t)| < 2e, n = 1,2, . . . . 
t>o 



The induction hypothesis is verified, using (438 1, (4.39 1, (4.42 1, and (4.44i: 



|a;2(i)| < eexp{— Kt} + 5e exp{— «;(t — s)}ds < 2e. 

Jo 



(4.48) 



Similarly, one shows that (4.48 )forn = k implies (4.48 ) for n = k + 1. Thus, (4.48 1 holds for all natural n. 



We complete the proof by showing that x„(t) uniformly converges the solution of (4.36 1 and (4.37 ). To 
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this end, we show 



sup|xfc+i(t) -Xkit)\ < -sup|xA;(t) - Xk-i{t)\, A; = 2,3, 



(4.49) 



Indeed, by subtracting ( |4.47[ ) with n = k from ( |4.47[ ) with n = k + 1 and using ( |4.38[ ), ( |4.45[ ), and ( |4.48| ), 
we have 



sup\xk+i{t) - Xk{t)\ < S \X(t,s)\sup\xk{t) - Xk-i{t)\ds < -su-p\xk{t) -Xk-i{t)\ 

t>0 Jo t>0 K t>0 



< -SUp\Xk{t) -Xk-l{t)\. 
^ t>0 



Next, consider 



^SUp|Xfc+l(t) - Xkit)\. 



(4.50) 



fc=i 



t>o 



By (4.49 1, (4.50 1 is majorized by the geometric series YlT=i ^ • Therefore, {x„(t)} converges uniformly 



to x{t), the unique solution of ( |4.36[ ), ( |4.37[ ). By ( |4.48| ), the latter (as the limit of {xn{t)}) is bounded by 2e. 
This completes the proof. 

D 



Proof. (Lemma 4.4 1 We continue to use the notation introduced in the proof of Lemma 4.5 In particular. 



positive e and 6 are as chosen above. By the variation of constants formula, we express the solution of ( 4.36 1 



and (4.37 1 as 



x{t) =X{t)xo+ X{t,s)Q{x{s))ds. 
Jo 



Using (4.38 1, (4.391, and (4.44i from (4.51 1 we have 



\x{t)\ < ex.p{—Kt}e + 6 / exp{—K{t — s)}\x{s)\ds. 

Jo 



(4.51) 



(4.52) 



Rewrite (4.52i for y{t) := \x{t)\ exp{«;i}: 



y{t) <e + 6 y{s)ds. 
Jo 
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By Gronwall's inequality, 



y{t) < eexp{5t}, 



and, by recalling the definition of y{t) and (4.42 ), we finally derive 



\x{t)\ < eexp{((5 — K)t} < eexp -^ — t 



D 
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Appendix. Parameter values for (3.9) and (3.10) 



The equations for the local systems in the neural network ( 3.9 1 and ( 3. 10 1 are adopted from a nondimensional 
model of a dopamine neuron |[36l (see also [341)- For biophysical background and details of nondimesion- 



alization, we refer an interested reader to II361 . Terms on the right hand side of the voltage equation (3.9 1 
model ionic currents: a calcium current, a calcium dependent potassium current, and a small leak current. 



The equation for calcium concentration (3. 10 1 takes into account calcium current and calcium efflux due to 
calcium pump. The ionic conductances are sigmoid functions of the voltage and calcium concentration 



gi{v) 



92[U) 



|fl + tanh^"-"^ 



02 



U^ + 4' 



Constants 51,2,3 and £'1,2,3 stand for maximal conductances and reversal potentials of the corresponding 
ionic currents; ai,2,3 are constants used in the descriptions of activation of calcium and calcium dependent 
potassium currents; lo and e are certain constants that come up in the process of nondimesionalization of the 
conductance based model of a dopamine neuron (see [36] for details). The values of parameters used in the 
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simulations shown in Figure |2] are summarized in the following table. 



Table 



El 


1 


E2 


-0.9 


Es 


-0.3 


51 


0.8 


92 


2 


93 


1 


9 


0.3 


ai 


-0.35 


a2 


1.4-10-2 


03 


1.8 


e 


0.1 


T 


5.0 


UJ 


5.0 
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